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Abstract. We study a contact process on a two-dimensional square lattice which is 
diluted by randomly removing bonds with probability p. For p < 1/2 and varying birth 
rate A the model was shown to exhibit a continuous phase transition which belongs to 
the universality class of strongly disordered directed percolation. The phase transition 
line terminates in a multicritical point at p = 1/2 and A = A* = 3.55(1), where the 
model can be interpreted as a critical directed percolation process running on a critical 
isotropic percolation cluster. In the present work we study the multicritical point and 
its neighboorhood by numerical simulations, discussing possible scaling forms which 
could describe the critical behavior at the transition. 
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1. Introduction 

Phase transitions from fluctuating phases into so-called absorbing states continue to be 
an important topic of non-equilibrium statistical physics. A prototypical system that 
exhibits such a non-equilibrium phase transition is the contact process. It is deflned 
on a c?-dimensional hypercubic lattice with sites that can be either empty (inactive) or 
occupied by a particle (active). As time evolves, particles are spontaneously removed 
at unit rate and created at empty sites at rate Xn/{2d), where the so-called birth rate 
A is a control parameter and n denotes the actual number of occupied nearest neighbor 
sites. In the inactive phase, where A is small, particle removal dominates and the 
system approaches a so-called absorbing state without particles from where it cannot 
escape. However, if A exceeds a certain threshold Ac, a steady state with a flnite density 
of particles exists on the inflnite lattice and the system is said to be in the active 
phase. The two phases are separated by a nonequilibrium phase transition at the critical 
threshold A = Ac- This transition is known to belong to the universality class of directed 
percolation (DP) [1-4]. 




Figure 1. Schematic pliase diagram of the diluted contact process as a function of 
the dilution parameter p and the birth rate A. The phase transition lines (a) and (b) 
are explained in the text. Both lines meet in the multicritical point (MCP) located at 
the percolation threshold p = pc and the critical birth rate A = A, . The present study 
focuses on the scaling behavior at and in the vicinity of the multicritical point. 



In spite of the progress in the theory of absorbing phase transitions, there is so 
far no experiment that reproduces the critical behavior of the transition in a rehable 
way [5] . Ahhough various experiments show signatures of a continuous transition on a 
quahtative level, none of them was able to reproduce the full set of DP critical exponents 
quantitatively. Designing and performing such an experiment remains one of the most 
challenging open problems of non-equilibrium statistical mechanics. 

There are various reasons why experiments of DP are so difficult to perform. For 
example, in the contact process the absorbing state is perfectly stable while in realistic 
situations there is always a small but finite rate for spontaneous creation of activity that 
washes out the transition. Another possible source of difficulties in many experiments 
is the presence of spatially quenched disorder due to inhomogeneities of the system. In 
the contact process this type of disorder can be introduced by initially varying the birth 
rate A from site to site and keeping these variations fixed as time evolves. From the field- 
theoretic point of view, spatially quenched disorder is a marginal perturbation under 
renormalization group transformations which entirely changes the critical behavior at 
the transition [6]. In fact, numerical studies by Moreira and Dickman [7,8] showed 
that the critical dynamics of the contact process with quenched disorder evolves on a 
logarithmic time scale, e.g., the density of particles decays as p{t) ~ [ln(t)]~'' with some 
exponent 6. More recently, Hooyberghs et al. [9, 10] argued that this type of critical 
behavior can be related to the random-field Ising model provided that the infiuence of 
quenched disorder is strong enough. 

A very simple model, where the effect of disorder can be studied systematically, 
is the diluted contact process introduced by Noest et al. [11, 12]. In this model a 
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two-dimensional square lattice is first diluted by randomly removing lattice sites with 
probability p. After diluting the lattice, a contact process runs on the remaining sites 
according to the usual dynamical rules. The schematic phase diagram of this model 
is shown in Fig. ^ For p = (no dilution) the model reduces to the ordinary two- 
dimensional contact process with a DP phase transition at Ac = 1.64877(3) [13]. For 
weak dilution below the percolation threshold of the lattice, p < Pc, the active phase of 
the model survives, but the critical birth rate A increases with p in order to compensate 
the missing neighbors. For strong dilution beyond the percolation threshold, p > pc, the 
lattice decomposes into many disconnected finite clusters so that a steady state with a 
positive density of active sites cannot exist. Right at the percolation threshold p = pc, 
the active phase survives only on the infinite percolation cluster provided that the birth 
rate A is larger than some threshold A*. Thus, the active phase is bounded by two 
different phase transition lines, a curved and a vertical one, denoted by (a) and (b) in 
the figure. The two lines meet in a multicritical point (MCP). 

Although the scaling properties along these transition lines (a) and (b) are already 
well understood, the scaling behavior at the multicritical point (MCP) has not been 
studied so far. At the multicritical point the model can be viewed as a critical DP 
process running on top of critical frozen clusters of the diluted lattice at the percolation 
threshold. As both models - isotropic percolation providing the support and the DP 
process running on top of it - are scale-free at the multicritical point and characterized 
by different critical properties, one expects a non-trivial interplay of both universality 
classes leading to an unconventional scaling behavior. The purpose of this paper is 
to present numerical results and to test possible scaling forms that could describe the 
critical behavior of the diluted contact process at and in the vicinity of the multicritical 
point. To this end the paper starts with a short review of various known types of scaling. 
In Sect.Elwe present our numerical results at and close to the multicritical point, testing 
several possible types of scaling. These results are then interpreted in the conclusions 
in Sect. H 

2. Conventional and non-conventional scaling forms 

In this section we briefly review several types of phenomenological scaling laws which are 
associated with different parts of the phase diagram in Fig. ^ A more comprehensive 
discussion can be found in Ref. [10]. 

2.1. Conventional scaling 

The clean contact process without dilution [p = 0) is known to exhibit ordinary power 
law scaling. This well-established type of scaling assumes that all relevant quantities 
such as distances, time intervals, the system size as well as order- and control parameters 
can be multiplied by a factor A raised to some power in such a way that the long-range 
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properties of the system close to the transition remain invariant. In the contact process 
the order parameter is the density of active sites p{A,t,L) which depends on the time 
parameter t, the distance from the critical point A = A — Ac, and the lateral system size 
L. In the vicinity of the clean transition (marked by 'DP' in the phase diagram), where 
A is small, this order parameter scales under changing lengths by a factor A as 

p(A, t, L) = A>(^A-i/^^ A, AH, Al) , (1) 

where x = P/i'± and z = vw/vi. are certain critical exponents. Scaling forms can be 
obtained by choosing A in such a way that one of the arguments on the r.h.s. becomes 
constant. For example, choosing A = one gets 

p(A,t,L)=t-V(t'/'^iiA,t-^/^L), (2) 

where 5 = x/z = (3/v\\ and / is a scaling function. In most cases this scaling function 
is universal, i.e., its functional form is fully given by the universality class to which the 
phase transition belongs. 

The above scaling form implies that the density of active sites measured in an 
infinite system at criticality decays as p{t) ~ t^^ . Moreover, approaching the transition 
the correlation length and the correlation time both diverge algebraically as 

U ~ A-^- , eii ~ A-^ii , (3) 

meaning that both scales are related by ~ 



2.2. Activated scaling 

Along the transition line (a) in Fig. ^ the diluted contact process is known to display 
an unconventional type of scaling behavior that has been termed as activated scaling 
or strong disorder scaling in the literature. This type of scaling is characterized by an 
ultra-slow temporal evolution which depends on In t instead of t. Activated scaling may 
be interpreted as a limiting case of conventional scaling when the exponent v\\ is taken 
to infinity. The resulting scaling laws are very similar to the conventional ones, the only 
difference being that t is replaced by Int and that the exponents may have different 
values. 

In the case of activated scaling the characteristic length scale still diverges 
algebraically whereas the correlation time diverges exponentially as the critical point 
is approached: 

a ~ A-^^ , eil - exp (aA-^) . (4) 

Here a is a constant and z7j[ is an exponent analogous to z/y in the conventional case. 
This means that time and length scales are no longer related by ^\\ = instead is 
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related to the logarithm of through 



(5) 



where z = The value of the exponents and z may depend on the strength of 

disorder, i.e., the site dilution probability p. 

In the case of activated scaling the usual relation for the density p upon rescaling 
L AL in Eq. ((H) is replaced by 



where x = (3/i>±_ and A denotes the distance from the transition line. Choosing 

A = [ln{t/to)]~^/^ one gets 



where F is a scaling function and 6 = x/l = (3 /V\\ plays a role of a decay exponent 
analogous to 5 for conventional scaling. Note that whenever one takes the logarithm of 
a parameter, it has to be divided by a constant to make the argument dimensionless. 
For example, in the scaling form given above, to plays the role of an elementary time 
scale. In a numerical simulation such a constant is expected to be of order 1 in natural 
units of Monte Carlo sweeps. 

As shown in a recent paper by Vojta and Lee [14], activated scaling is not only 
valid along line (a) but also along the vertical transition line denoted by (b) in FigQ 
Here the supporting lattice is at the percolation threshold and thus decomposes into 
many clusters according to a scale-free distribution. Each of these clusters hosts an 
independent contact process. In the limit t — > cxo a nonzero steady-state density can 
only survive on the infinite percolation cluster while finite clusters enter the absorbing 
state by rare fluctuations and thus do not contribute asymptotically. 

Since the infinite cluster has a fractal dimension Df < 2 and therefore covers only 
a vanishing fraction of lattice sites in an infinite system, it does not contribute to the 
density of active sites p(t). This means that p(t) is dominated by contact processes on 
finite clusters and hence approaches zero as t goes to infinity. However, this decay is 
extremely slow since contact processes on larger clusters may have a very long life time. 
Vojta and Lee showed that the transition along line (b) obeys activated scaling, i.e., we 
can use the same scaling form as for strongly disordered DP, although with a different 
set of exponents. 

Vojta and Lee also proved that the exponents (3 and z/^ on line (b) coincide with 
the usual order parameter and correlation length exponents jSc = 5/36 and z/^ = 4/3 
of isotropic percolation in two dimensions, i.e., all static features are dictated by the 

I In previous works the exponent z is also denoted as ip. 




(6) 
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scaling properties of the supporting lattice at the percolation threshold. For example, 
in the vicinity of the critical line the steady-state density p^*"* and the correlation length 
are scale as 

P^*"*~b-Pci''% u-\p-Pcn- {p<Pc) (8) 

Moreover, they found that the density of active sites decays as 

p(t) ~ [ln(^/^or^ (9) 

where 6 = Pc/{^cDf) = 5/91 is again related to the exactly known static exponents and 
the fractal dimension of isotropic percolation in two dimensions. In the vicinity of the 
line (b) sufficiently far away from the multicritical point one finds a Griffiths-like power 
law decay 

Pit) ~ (t/to)-"/'' (10) 
for p > pc and a stretched exponential of the form 



pit) - pst ~ exp 



d rt^^'-'/' 

m 

to 



:iii 



for p < Pc, where pst again denotes the stationary density. Here d is the spatial dimension 
while z' and z" are nonuniversal exponents with diverge as p —>■ Pc- 

To summarize, both transition lines (a) and (b) obey activated scaling, although 
with different sets of critical exponents, the latter being related to the exactly known 
exponents of isotropic percolation in two dimensions. However, the physical mechanisms 
leading to a logarithmic time dependence are very different in both cases. As pointed 
out by Vojta and Lee, this difference is also reflected in the early time behavior of the 
model. To this end one considers the dynamics of a contact process starting from a 
single active site on a diluted lattice. On line (b), where the contact process is locally 
supercritical, the cloud of active sites first grows linearly with time until it covers its 
supporting cluster completely. This cloud then remains in a metastable state until it 
reaches the absorbing state by a rare fiuctuation. Therefore one observes a very fast 
short-time behavior followed by a logarithmically slow decay. On line (a), however, the 
observed initial spreading is much slower. 



2.3. Scaling in the Griffiths phase 

Quenched disorder by dilution leads to the emergence of a so-called Griffiths phase which 
in the present model is located between the transition line and the horizontal line starting 
at the clean critical point (the dashed line in Fig^. A Griffiths phase is characterized 
by an algebraic decay of the order parameter with non-universal continuously varying 
exponents while spatial correlations remain short-ranged. 
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The algebraic decay can be understood as follows. Consider a compact finite region 
of the lattice and let k be the corresponding number of bonds. The probability that 
none of these bonds was removed during the dilution procedure is (1 —p)'', hence the 
probability to find a fully connected region decays exponentially with its size. Since for 
large regions with typically more than — ln(l — p) bonds this probability is extremely 
small, such fully connected islands have been named as rare regions in the literature. 
For A > Ac each rare region supports a locally supercritical contact process that survives 
in a metastable active state for a typical time r which grows as r oc exp(+Q;A;) with 
k, where a > depends on A — Ac. The exponentially diverging life time combined 
with an exponentially decaying likelihood for rare regions is the origin of the observed 
non-universal power-law behavior in the Griffiths phase. 

To see this note that at a given time t only those rare regions are still in the 

metastable active state whose number of bonds k is typically larger than a Int. Summing 
up their contributions to the particle density weighted by the probability (1 — p)*^ to 
find such regions, one arrives at an expression proportional to 

(l_p)(lnt)/a ^ ^ln(l-p)/a_ ^^2) 

Obviously, the continuously varying exponent ln(l — p)/q; depends on both p and A — Ac. 



2.4- Scaling forms based on continuously varying exponents 

In the present work we also test another type of logarithmic scaling which has been 
discussed in the context of multiscaling. The theoretical background for this generalized 
scahng theory is described in Ref. [15] and can be sketched as follows. 

Starting point is to generalize conventional scale transformations 

T(a) : {A ^ A-V-x A , t ^ AH , L ^ AL , A'^/^^p} (13) 

in such a way that the critical exponents (3, z/_|_, z/y are replaced by continuously varying 
exponents which may depend on the parameters A, t and L. However, the functional 
form of these exponents is not arbitrary, rather it is restricted by the postulate that 
T(a) can still be interpreted as a scale transformation (dilatation) controlled by a scale 
parameter A. More specifically, it is required that two subsequent scale transformations 
by Ai and A2 are equivalent to a single transformation by the factor A1A2, i.e., scale 
transformations should obey the group homomorphism 

t(Ai) o r(A2) = t(AiA2) . (14) 

As shown in Ref. [15], this requirement leads to the partial differential equations 



D[P{A, t, L)] = D[u^{A, t, L)] = D[u\\ (A, t, L)] = 



(15) 



Multicritical behavior of the diluted contact process 



8 



with the differential operator 

O O r\ 

D = P{A,t,L)A— + u\\{A,t,L)t- + u^{A,t,L)L—. (16) 

These differential equations constrain the possible functional form of the continuously 
varying critical exponents. Since D is the generator of generalized scale transformations, 
these partial differential equations tell us that the exponents may vary with the 
parameters but only in such a way that they are invariant under scale transformations r. 
Clearly, conventional scaling with constant critical exponents is a possible solution. 
However, as shown in Ref. [15], there are also other non-trivial solutions. For simplicity 
let us consider the critical case A = 0, where p(t, L) is a function of only two parameters. 
One possible solution assumes that the exponents and the density are functions of the 
quotient 

ln(t/to) 
ln(L/Lo) 

and that the exponents z^||(a;) and v±{a) are related by 

z/||(a) = az/_L(a). (18) 
In order to verify this solution, note that Eq. p8p turns the generator D into 

SO that Dh{a) = for any function h that depends exclusively on a. 

In addition, the functional form of the order parameter p{t, L) is restricted by the 
differential equation Dp(t,L) = P{a)/h'j_{a). Solving these differential equations one 
arrives at the scaling form 

which can be tested numerically by an appropriate data collapse. Again we divided all 
relevant quantities by constants to make the arguments of the logarithm dimensionless. 
Setting H = this logarithmic scaling form was successfully apphed to several systems 
that are believed to exhibit multiscaling, including self-organized critical sandpile 
models [16,17], DLA-related growth processes [18,19] and experiments in turbulence [20]. 
However, in spite of its apparent generality the logarithmic scaling form seems to fail in 
the present will be shown. 
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t=l t=10 t=100 t=10000 



Figure 2. Snapshots of the bond-diluted contact process on a 300 x 300 lattice with 
periodic boundary conditions at the multicritical point. Active sites are represented 
by red pixels. The maximal cluster, which in a numerical simulation approximates the 
infinite one, is highlighted by green pixels. 



3. Numerical results 

3.1. Definition of the model 

The diluted contact process studied in this paper is defined on a two-dimensional square 
lattice with periodic boundary conditions and N = LxL sites which can be either active 
or inactive. In contrast to Refs. [11,12,21], we will consider a bond-diluted contact 
process, i.e., bonds instead of sites are randomly removed with probability p. This type 
of dilution has the advantage that the percolation threshold of the lattice pc = 1/2 is 
known exactly [22]. 

After removing the bonds, all lattice sites are initially occupied by a particle. The 
model then evolves by random-sequential updates in the same way as the ordinary 
contact process: For each attempted update one of the particles is randomly selected. 
With probability this particle is removed, otherwise one of the four nearest 
neighboring sites is randomly selected and a new particle is created provided that the 
target site is empty and the corresponding bond has not been removed during the initial 
dilution procedure. Each update attempt corresponds to an average time increment of 
' where n is the actual number of active sites. Typical snapshots of a simulation of 
a single run are shown in Fig. [2l The quantities to be measured are averaged over a large 
number of runs, each of them with an independent realization of quenched randomness. 

As usual, the simplest order parameter is the average density of active sites 
p{t) = n/N. Without dilution [p = 0) at the DP critical point A = Ac, this order 
parameter is known to decay algebraically as p{t) ~ t~^ on an infinite lattice, where 
5 = ~ 0.4505 is the corresponding critical exponent. In what follows we study the 
diluted contact process along the line p = 1/2, where the frozen network of bonds is at 
the percolation threshold. 
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Figure 3. Stationary density of active sites on the maximal cluster at the percolation 
threshold p ^ Pc for various A. The solid line shows a fitted power law (see text). 

3.2. Steady-state properties 

The static exponent /3 can be estimated by measuring the stationary density of active 
sites close to the transition. However, keeping p = 1/2 fixed and varying A one finds 
that the stationary density of active sites p**"* is always zero, even on line (b), where the 
contact process is locally supercritical. This is because the only surviving subprocess at 
the percolation threshold is the one running on the infinite percolation cluster while all 
other subprocesses on finite clusters have a finite (although sometimes very long) life 
time. As the infinite percolation cluster with Df < 2 covers only a vanishing fraction of 
the lattice, this surviving subprocess does not contribute to the overall density of active 
sites. 

In numerical simulations, however, one finds a quite stable positive density of active 
sites. The reason is that in finite system the maximal cluster still covers a considerable 
finite fraction of the lattice and therefore contributes to the overall density of active 
sites. As the fractal dimension Df = 91/48 ~ 1.896 of critical isotropic percolation is 
close to 2 these finite-size corrections are enormous, e.g. on a lattice of 2000 x 2000 sites 
the largest cluster typically occupies about 30% of the lattice. 

To circumvent this problem, we study the density of active sites on the maximal 
cluster (which in a simulation approximates the infinite one) . To this end we first dilute 
the lattice by removing bonds with probability p, classify all clusters using of the Hoshen- 
Kopelman algorithm [23], identify the largest cluster and monitor the contact process 
exclusively on this cluster during the subsequent temporal evolution. The corresponding 
density, defined as the number of active sites on this cluster divided by its mass, will be 
denoted as pmi'^,t). 
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Figure 4. Decay of the density of active sites on the maximal cluster pm{t) at the 
multicritical point p — Pc = 1/2 and A = A* = 0.355 versus t and In t in a double- 
logarithmic representation. The dashed line indicates the slope -1.63. 



Fig. |H1 shows the stationary density on the maximal cluster function of A 

measured in a system of size 2000 x 2000 averaged over 100 runs. Although the errors 
of the estimates are quite large, especially close to the critical point, the results are 
compatible with power-law scaling of the form 

Pm^Ai\-K)^, (21) 

with 

A,, = 3.55(1), /3 = 0.81(7) (22) 
and the proportionality factor A ^ 0.464. 



3.3. Temporal decay at the multicritical point 

Having determined A^. we study the temporal decay of Pmit) at the multicritical point. 
The numerically determined curve is plotted in Fig. |3] in two different representations. 
The left panel shows the data in an ordinary double- logarithmic plot. As can be seen, 
the line has an overall curvature, hence conventional power-law scaling is very unlikely. 
However, if the same curve is plotted against ln(t) in a double logarithmic plot, one 
obtains after an initial transient a straight line extending over almost four decades in 
time, as shown in the right panel of Fig. HI This suggests a logarithmically slow decay 
of the form 

Pm{t) ~ Mt/toT' (23) 

with an exponent 6 = 1.63(10) and the constant to = 1- We note that the choice of 
to is crucial for the estimate of the exponent (see discussion in Sect. H}. In the present 
analysis we used to = 1 for which the cleanest logarithmic decay is observed. 
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Figure 5. Off-critical simulations close to the multicritical point at the percolation 
threshold p — 1/2. Left: Temporal decay of the density of active sites measured on the 
maximal cluster in a system with 2000 x 2000 sites for various values of the birth rate 
A. The curve for X — X^, where the density decays logarithmically with time, is marked 
by a bold line. Right: Data collapse according to the scaling form (|24(l for A varying 
from 3.3 to 3.9 in steps of 0.05. In order to avoid non-universal initial transients, this 
collapse shows only data points after 10 Monte-Carlo sweeps. 



3.4- Off-critical scaling 

Next we study the temporal decay of Pm{^, t) for A = A — A* 7^ in the vicinity of the 
multicritical point, keeping p = 1/2 fixed. Our numerical results are shown in the left 
panel of Fig. El As can be seen, for A < A* (A < 0) one has straight lines with varying 
slopes, i.e. the system crosses over to an algebraic decay with continuously varying 
exponents which is the typical behavior in a Griffiths phase. For A > A* (A > 0), i.e. on 
line (b) of the phase diagram, the contact process on the maximal cluster is supercritical 
so that the curves saturate at a constant value p^'^*(A). 

Since the stationary density was shown to increase with A as a power law while at 
criticality the decay is logarithmically slow, it is very likely that the order parameter 
Pm{^, t) exhibits activated scaling. The corresponding scaling form can be derived from 
Eq. (jni) by inserting A = A"^ and taking L to infinity, giving 

p™(A,t) = A'^F(A^ln(t/to)). (24) 

Hence, plotting pm{^, t)/Sr^ versus A"^!! ln(t/to) this scaling form can be tested by a data 
collapse. As shown in the right panel of Fig. the supercritical curves along line (b) 
can be collapsed convincingly for 17\\ = 0.50(5). The subcritical curves in the Griffiths 
phase, however, do not collapse entirely. This may be related to the non-universality of 
the power-law decay in the Griffiths phase. 
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Figure 6. Decay of the density of active sites for (a) an contact process on a finite 
lattice with TV = 32^ = 1024 sites at A = A™ averaged over 27000 runs compared to 
(b) a contact process on a percolation cluster with approximately 1024 sites at A = A,, 
averaging over 35000 independent clusters. In both cases the density p{t) is multiplied 
by t'. 



3.5. Finite-size scaling 

So far we analyzed the scaling properties of the density of active sites pm(A,t) on the 
maximal cluster, which in a numerical simulation can be regarded as approximating the 
infinite cluster at the percolation threshold. Let us now turn to the remaining clusters. 
Since in two-dimensional critical percolation the cluster sizes are distributed as [22] 



most of these clusters are very small so that finite-size effects are expected to play an 
important role. 

As illustrated in Fig. finite-size effects on finite percolation clusters differ 
significantly from ordinary finite-size effects. As an example the figure shows the decay 
of the density of active sites for (a) a critical contact process on a finite lattice with 
32^ = 1024 sites, compared to (b) a contact process running on a percolation cluster 
with approximately 1024 sites, averaged over many independent runs. As can be seen, 
the resulting curves differ significantly. Qualitatively this difference can be explained as 
follows: 

(a) On a finite lattice the system first evolves in the same way as in an infinite system, 
exhibiting a clean power-law decay p(t) ~ t^^ and a growing correlation length 
soon as C,±{t) becomes comparable with the lateral size L of the 
system (32 sites in Fig. the initial power law crosses over to an exponential 
decay. As can be seen in the figure, this finite-size effect sets in quite suddenly at 
a well-defined typical time. 




(25) 
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[ln(t/tQ)] / s 

Figure 7. Finite size scaling analysis. Left: Density active sites for various cluster 
sizes s ranging from 5 to 4 x 10® and grouped geometrically into 10 bins per decade. 
Right: Data collapse according to the scaling form H26|) . as explained in the text. 



(b) On a percolation cluster with about the same number of sites, finite-size effects 
seem to influence the decay of the particle density at any time t, leading to a soft 
overall curvature of the density in the double- logarithmic plot. This is plausible 
since on a percolation cluster there is no well-defined distance from boundary to 
boundary, instead the fringes of the cluster induce a whole spectrum of length 
scales which have to be compared with the growing correlation length of the contact 
process. Therefore, finite-size effects are always present and become gradually more 
important as time proceeds. 

To analyze the decay of the particle density on finite clusters at the multicritical point, 
we first classify all clusters at the beginning of each run using the Hoshen-Kopelman 
algorithm [23]. For the subsequent statistical analysis the cluster sizes s (i.e., the number 
of sites occupied by the clusters) are grouped logarithmically into bins with a density 
of 10 intervals per decade. For our simulations on a lattice with 2000 x 2000 sites this 
requires about 67 bins. For each bin we monitor the activity of all corresponding clusters 
as a function of time, averaging over 100 independent runs. In this way we obtain a set 
of 67 curves parameterized by s, each of them describing the temporal evolution of the 
density of active sites on clusters with typical size s. 

In analogy to the off-critical case we test the possibility of activated scaling for finite 
clusters by a data collapse. Here the lateral system size L in Eq. (jH)) has to be replaced 
by the typical lateral length scale ^ of a finite percolation cluster in two dimensions. For 
a cluster with s sites this length scale is expected to be proportional to £ ~ s^^^f , where 
Df = 91/48 ~ 1.896 is the fractal dimension of isotropic percolation in two dimensions. 
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Inserting A = i ^ = s ^^^f into Eq. (jH)) one obtains the scaling form 

p{s, t) = s-^^/^f F (^s-'/^f ln(t/to)) • (26) 

Hence, plotting p(s,t)s^^^'^^^f^ versus ln(t/to)/s^^^^ this scaling form can be tested 
numerically by a data collapse. As shown in Fig.[7|a convincing data collapse is obtained 
for a = 'z/Df = 0.30(2), hence we arrive at the estimate 

z = 0.57(4) . (27) 
3. 6. Failure of logarithmic scaling 

Finally we demonstrate that generalized scaling with continuously varying exponents 
according to Ref. [15] seems to fail in the present model. As outline in Sect. 12. 4[ this 
theory predicts a logarithmic scaling form which in the present case reads: 



^ ^/ ln(A/Ao) ^ ln(tAo) \ 
\ ln(s/so) ' ln(s/so)/ 



Unlike activated scaling, this type of scaling requires four normalization constants to, -Sq 
and po which may be regarded as fit parameters. Typically one expects these constants 
to be of the order 1 in natural units of the simulation. However, even when these 
constants are varied, we do not find convincing data collapses for both off-critical and 
finite-size simulations. For H = the failure of the data collapses is demonstrated in 
Fig. IHl where we set Aq = to = sq = po = 1. 

There is also a theoretical hint why this scaling form fails: Taking A — > and 
s ^ oo one can show that Eq. (^Hj) reduces to l^^tjl^^^ = const, i.e., p{t) has to decay 
algebraically provided that this constant is nonzero. As demonstrated in Sect. 13. 3| this 
is not the case, instead a logarithmic decay is found. 
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transition 


quantity 


scaling type 


/3 




z,z 


S,S 


Ref. 


DP (2d) 


P, Pm 


conventional 


0.556(5) 


0.76(1) 


1.766(1) 


0.451(1) 


[4] 


line (a) 


P, Pm 


activated 


0.93(5) 


1.00(9) 


1.98(12) 


0.47 


[8] 


line (b) 


P 


activated 


5/36 


4/3 


91/48 


5/91 


[14] 


MCP 


Prn 


activated 


0.81(7) 


0.88(10) 


0.57(4) 


1.63(10) 


this work 



Table 1. Numerical estimates of the critical exponents. 



4. Discussion 

In the present paper we have analyzed the scahng behavior of the bond-diluted contact 
process at the percolation threshold p = 1/2 for various birth rates A in the vicinity of the 
multicritical point A = A* = 0.355(1). Along this line the diluted lattice decomposes into 
a large number of clusters according to a scale-free distribution, each of them hosting 
an independent contact process. Performing numerical simulations we arrive at the 
following results: 

(i) The overall density p, which monitors the average activity of all contact processes, is 
characterized by a complex mixture of finite-size effects. While the infinite cluster 
does not contribute to p in an infinite system, it leads to severe corrections in 
numerical simulations even even if one uses very large lattice sizes. Therefore, as a 
more appropriate order parameter, we studied the density of active sites restricted 
to the maximal cluster, Pm{t), which is non-zero along line (b) and thus allows us 
to define a static exponent (3. 

(ii) The temporal decay of this order parameter at the multicritical point Pm{t) ~ 
[ln(t/to)]^'^ is logarithmically slow, suggesting activated scaling at and in the vicinity 
of the multicritical point. This conjecture is supported by off-critical simulations 
and by an analysis of the density decay on finite-size clusters at criticality. 

(iii) The estimated exponents (see Table differ from those on line (a) and line (b), 
suggesting a crossover phenomenon similar those observed at ordinary multicritical 
points with power-law scaling. 

(iv) A fully logarithmic scaling form suggested in Ref. [15] seems to fail in the present 
model. 

As a cross-check we may now use these results in order to predict the decay of the 
overall density p{t) at the multicritical point. To this end one has to integrate p{s, t) over 
all cluster sizes, weighted by the cluster size s and the probability P{s) ~ s^^ = s^i87/8i 
to find such clusters: 

p{t) = J dssP{s) [ln(t/to)]-'F([ln(t/to)]/s") • (29) 

Dimensional analysis of this integral yields 



Pit) ~ [ln(t/to)]-^-(^'^)/' 



(30) 
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Figure 9. Overall density of active sites p{t) compared to the density on the maximal 
cluster pm{t) versus t. As can be seen, p{t) decays faster then Pm{t)- The inset shows 
the same curves versus In i in a double- logarithmic representation. The slopes of the 
dashed lines indicate the theoretical prediction according to Eqs. H23|l and (|30|l . 

This means that the overall density decays faster than the density on the maximal 
cluster. In fact, this prediction can be verified numerically, as demonstrated in Fig. IHl 

Let us conclude with a critical remark. Unlike previous studies along lines (a) 
and (b), the results of the present work are based entirely on numerical evidence. 
Although the hypothesis of activated scaling at the multicritical point seems to be 
plausible, we realized that numerical simulations with an ultra-slow temporal evolution 
on a logarithmic time scale are very difficult to perform. In particular they are much 
more susceptible to misinterpretations caused by a wrong expectation. For example, 
if we expected power-law behavior at the transition, we would find an almost perfect- 
looking straight line in a log-log plot at p = 1/2 and A = 3.47. Even data collapses 
in off-critical and finite-size simulations could be generated, but compared to the ones 
presented above their quality would be 'less convincing'. The suggested scenario of 
activated scaling is much more coherent but still far from a proof. 

Another important source of incertitude in scaling laws involving logarithms 
concerns the normalization constants needed to make arguments dimensionless. These 
constants are expected to be non-universal and have to be treated as fit parameters. For 
activated scaling the crucial constant is the time scale to- In our analysis we set t^ = 1 
for which the cleanest logarithmic decay of the density is seen. But obviously there is no 
particular reason for this particular choice, e.g. a simple redefinition of the rates in the 
Monte Carlo procedure would require a corresponding change of to- We were unable to 
specify a reliable error margin for to- For example, setting to = 1/2 the data collapses 



Multicritical behavior of the diluted contact process 



18 



would be slightly less convincing but the whole analysis still works but the estimates of 
the critical exponents seem to depend on the choice of t^. This source of errors is not 
yet satisfactorily controlled. 
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